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Well-balanced finite volume evolution Galerkin 
methods for the shallow water equations 


M. Lukacova - Medvid’ov^, S. Noell^ and M. Kraft^ 


Abstract 

We present a new well-balanced finite volume method within the framework of the 
finite volume evolution Galerkin (FVEG) schemes. The methodology will be il¬ 
lustrated for the shallow water equations with source terms modelling the bottom 
topography and Coriolis forces. Results can be generalized to more complex sys¬ 
tems of balance laws. The FVEG methods couple a finite volume formulation with 
approximate evolution operators. The latter are constructed using the bicharacter¬ 
istics of multidimensional hyperbolic systems, such that all of the infinitely many 
directions of wave propagation are taken into account explicitly. We derive a well- 
balanced approximation of the integral equations and prove that the EVEG scheme 
is well-balanced for the stationary steady states as well as for the steady jets in the 
rotational frame. Several numerical experiments for stationary and quasi-stationary 
states as well as for steady jets confirm the reliability of the well-balanced EVEG 
scheme. 


Key words: well-balanced schemes, steady states, systems of hyperbolic balance laws, 
shallow water equations, geostrophic balance, evolution Galerkin schemes 
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1 Introduction 

Consider the balance law in two space dimensions 

Ut + fiiu)^ + f 2 iu)y = b{u, X, y), (1.1) 

where u is the vector of conservative variables, /l, /2 are flux functions and b{u,x,y) is 
a source term. In this paper we are concerned with the hnite volume evolution Galerkin 
(FVEG) method of Lukacova, Morton and Warnecke, cf. [22]-[29] and [I5]. The FVEG 
methods couple a hnite volume formulation with approximate evolution operators which 
are based on the theory of bicharacteristics for the hrst order systems [23]. As a result 
exact integral representations for solutions of linear or linearized hyperbolic conservation 
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laws can be derived, which take into account all of the inhnitely many directions of wave 
propagation. 

In the finite volume framework the approximate evolution operators are used to evolve 
the solution along the cell interfaces up to an intermediate time level tn+i /2 in order to 
compute fluxes. This step can be considered as a predictor step. In the corrector step 
the hnite volume update is done. The FVEG schemes have been studied theoretically as 
well as experimentally with respect to their stability and accuracy. Extensive numerical 
experiments confirm robustness, good multidimensional behaviour, high accuracy, stabil¬ 
ity, and efficiency of the FVEG schemes, see Section |5] and also references [26l [27] . We 
refer the reader to [iiEiiTiiHiigi usi nniiao] for other recent multidimensional schemes. 

For balance laws with source terms, the simplest approach is to use the operator splitting 
method which alternates between the homogeneous conservation laws 

Ut + fi{u)^ + f2{u)y = 0 
and the ordinary differential equation 


Ut = b{u,x,y) 


at each time step. For many situations this would be effective and successful. However, the 
original problem fll.ip has an interesting structure, which is due to the interplay between 
the differential terms and the right-hand-side source term during the time evolution. For 
many flows which are of interest in geophysics, the terms are nearly perfect balanced. 
If these terms are treated separately in a numerical algorithm, the fundamental balance 
may be destroyed, resulting in spurious oscillations. In particular, we will be interested 
in approximating correctly equilibrium states or steady states, for which 

fl{u)x + f2{u)y = b{u,x,y), 

and we want to approximate perturbations of such equilibrium states. Equilibrium solu¬ 
tions play an important role because they are obtained usually as a limit when time tends 
to infinity. 

In this paper we present an approach which allows to incorporate treatment of the source 
in the framework of the FVEG schemes without using the operator splitting approach. 
The key ingredient is a new approximate representation of the multi-dimensional solution 
which contains the full balance of the hydrostatic pressure and the source terms, see 
Lemma 13.21 This new representation allows to apply a recent non-standard quadrature 
rule (see fl3.21l) and also [27|) to the mantle integrals of the bicharacteristic cone. This 
leads to a surprisingly simple, accurate and efficient approximate evolution operator and 
hence to an efficient, accurate and well-balanced finite volume scheme. We call our scheme 
the well-balanced finite volume evolution Galerkin scheme (FVEG). We refer the reader 
to OmEl UniHSlimlini eh uni IM] and the references therein for other related approaches 
for well-balanced finite volume and hnite difference schemes. 

Our paper is organized as follows. In Section 2 we introduce the shallow water equations 
in a rotating frame, derive a class of equilibria which contains the lake at rest, jets in 
the rotating frame and combinations of these two solutions. Then we introduce the 
class of two-step hnite volume schemes used throughout the paper. In Theorem 12.11 
we give sufficient conditions which guarantee well-balancing. Section [31 is devoted to 







Well-balanced EG schemes. 


Revised version, May 31, 2006 


3 


the EG (evolution Galerkin) predictor step. Applying the theory of bicharacteristics for 
multidimensional hrst order systems of hyperbolic type the exact evolution operator is 
derived. A well-balanced approximation of the exact evolution operator, which preserves 
some interesting steady states exactly and also works well for their perturbations, is given 
in Lemma and equations fl3.25p - fl3.26p . In Theorem 13.11 we prove that the FVEG 
scheme is well-balanced for the stationary steady states (e.g. lake at rest), for the steady 
jets on the rotating plane as well as for combinations of these two flows. In Section 4 we 
summarize the main steps of the FVEG method and present its algorithm. Numerical 
experiments for one and two-dimensional stationary and quasi-stationary problems as well 
as for steady jets presented in Section 5 conhrm reliability of the well-balanced FVEG 
scheme. We also compare the accuracy and runtime of the new FVEG scheme with that 
of the recent high order well-balanced WENO FV schemes of Audusse et a. [2] and Noelle 
et al. [31]. The question of positivity preserving property of the scheme, i.e. h > 0, is not 
yet considered here and will be addressed in our future paper. 


2 Geophysical equilibria and well-balanced two-step 
finite volume schemes 


2.1 The shallow water equations 


There are many practical applications where the balance laws and the correct approxi¬ 
mation of their quasi-steady states are necessary. Some example include shallow water 
equations with the source term modelling the bottom topography, which arise in oceano¬ 
graphy and atmospheric science, gas dynamic equations with geometrical source terms, 
e.g. a duct with variable cross-section, or fluid dynamics with gravitational terms. In 
what follows we illustrate the methodology on the example of the shallow water equations 
with the source terms modelling the bottom topography and the Goriolis forces. This 
system reads 


Ut /i('u)x + /2('“)y = 


where 


u 


AN 


h \ 

hu , /i(u) 

hv J 


hu 

hu^ + ^gh'^ 
huv 




hv \ 
huv , b{u) 

hv^ + \gh^ ) 


0 \ 

-ghbx fhv 

—ghby — fhu j 


( 2 . 1 ) 


Here h denotes the water depth, u, v are vertically averaged velocity components in x— 
and y— direction, g stands for the gravitational constant, / is the Goriolis parameter, and 
6(a;, y) denotes the bottom topography. 

Note that these equations are also used in climate modelling and meteorology for geostrophic 
flow, see, e.g., mm- For simulation of river or oceanographic flows some additional terms 
modelling the bottom friction need to be considered as well. 
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2.2 Equilibria 

Many geophysical flows are close to eqnilibrium, or stationary state. It is easiest to identify 
these states when the system is written in primitive variables. 


Wt + Ai{w)w^ + A2{w)wy = t{w), 


h \ / u h 0 \ / V 0 h 

w = \ u j,Ai=j g u 0 j,A2=j 0 v 0 
V j yOOny \ g 0 V 

Here we consider states which are both stationary, 

{h,u,v)t = 0 , 


(^)l 

(2.2) 

( ° \ 


, i = -gbx + fv 

. (2.3) 


-gby - fu 


(2.4) 


and constant along streamlines, 

(h, ft, i}) = 0, (2.5) 

where ' = dt + udx + vdy is the material derivative. For snch states, we obtain 

uhx + vhy = uux + vuy = uvx + vvy = 0 . ( 2 . 6 ) 

Since the velocity vector is now constant along the streamlines, these become straight lines. 
It is then natural to align the coordinates with the streamlines. The desired solution has 
to satisfy the conditions 


u 

= 0 

(2.7) 

Vy 

= 0 

(2.8) 

vhy 

= 0 

(2.9) 

g{h -h b)x 

= fv 

(2.10) 

g{h + b)y 

= 0. 

(2.11) 


In the region {{x,y) \ v(x) = 0} we obtain the lake at rest solution, where the water level 
h + b is flat. When v(x) ^ 0, we must have hy = Q and hence by = 0. Hence the 
topography is locally one-dimensional, and along the rise of the bottom we have a one¬ 
dimensional flow. This solution, which is well-known to oceanographers, is called the jet 
in the rotational frame. Due to the earths rotation the jet exerts a sidewards pressure fv 
onto the water, which is balanced by a raise in the water level g{h + b)x. In meteorological 
literature this state is also called the geostrophic equilibrium. 

For future reference we also dehne the primitive {U, V) of the Coriolis force, as introduced 
by Bouchut et al. in [2], via 


f f 

Vx = —V and Uy = —u 
9 " 17 


and the potential energies 


( 2 . 12 ) 


K\=g{h + b — V) and L ■.= g{h + b + U). 


(2.13) 
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2.3 Two-step finite volume schemes 

The FVEG (Finite Volume Evolution Galerkin) scheme which we propose for the balance 
laws are time-explicit two-step schemes, similarly as Richtmyer’s two-step version of the 
Lax-Wendroff scheme [32] and Golella’s GTU (corner transport upwind) scheme |3]. 

The hrst step, called predictor step, evolves the point value at a quadrature node to the 
half-timestep. This can be done by a simple hnite difference operator as in [32], by one¬ 
dimensional characteristic theory as in [6] or by fully multidimensional, bicharacteristic 
theory as in [25] and related works of Lukacova, Morton, Warnecke et ah. Our predictor 
step is based on this bicharacteristic theory. 

The second step is the standard hnite volume update. It approximates the hux integral 
across the interfaces by a quadrature of the huxes evaluated at the predicted states at the 
half-timestep. 

We proceed as follows: in the present section we study the hnite volume step (i.e. the 
second of the two steps). This will give us sufficient conditions for well-balancing which 
should be satished by the values computed in the predictor step (the hrst step). In 
Section 13.11 we will introduce the evolution Galerkin predictor step and prove that it 
satishes the sufficient conditions derived in the present section. 

In order to dehne the class of two step hnite volume schemes, let us divide a computational 
domain G into a hnite number of regular hnite volumes flij = [Xj_i, x , 2/j+i] = 
[xi — h/2, Xi + h/2] X [uj — h/2, Hj + h/2], i,j G Z, where h is the mesh size. Denote by 
the piecewise constant approximate solution on a mesh cell Qij at time tn and start with 
initial approximations obtained by the integral averages C/^- = U(-,0). Integrating 

the balance law (12.11) and applying the Gauss theorem on any mesh cell Hij yields the 
following hnite volume update formula 


u: 


n-l-1 


= u 


V 




n-l-l/2 


+ XB 


n-l-l/2 




k=l 


(2.14) 


where A = At/h, At is a time step, 6'^J^ stands for the central diherence operator in the 
Xfc-direction, k = 1,2 and represents an approximation to the edge hux at the 

intermediate time level tn + At/2. Further stands for the approximation of the 

source term multiplied with the mesh size, hb. The cell interface huxes are evolved 

using an approximate evolution operator denoted by EAi /2 to f„-|-At/2 and averaged along 
the cell interface edge denoted by £, 


^..+1/2 (2,15) 

j 

Here x^{£) are the nodes and Uj the weights of the quadrature for the hux integration 
along the edges. 

For simplicity, we introduce the following notation. Along the edges, we have quadrature 
nodes (a^j±i,2/j+y) resp. {xi+ti where i',j' G {0,±|}. These nodes are already 
sufficient for the midpoint, the trapezoidal and Simpson’s rule. We denote the values 
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given at the predictor step by 



j+i' 



and 


j+j' 





and the corresponding fluxes in x resp. ^/-direction by 

f ■ f 1 

^2 

With this notation, we obtain 




rij ^^^+1/2 
^X2J 2 


j' 







Finally we discretize the source term by 

/ 


^n+ 


-Bj, " = - J 


\ 


Ej. ikj+i') (b - r)«+/) 

Ef‘*'1' fci+i'j) (iSit'’' (b + B)i+i'j) / 


(2.16) 

(2.17) 

(2.18) 

(2.19) 

( 2 . 20 ) 


( 2 . 21 ) 


Here U and V are the discrete primitives of the Coriolis forces (see fl2.12|) l. dehned by 


y*0+i'T>. . — kL 

X*+*'j77. ,, . _ fcZ //*+*'’R'',. . 

^X 2 ^ h^X 2 


and we have used the average operators 


h-xi® ~ (®j+l/2,j + flj-l/2,j)/2 
“ (®ij+l/2 + flij-l/2)/2. 


( 2 . 22 ) 

(2.23) 


The following theorem states conditions which guarantee that the two-step finite volume 
scheme (I2.14p ~ (12.151) is well-balanced for the lake at rest as well as for the jet in the 
rotating frame. 

Theorem 2.1. Suppose that the values {h,u,v) given by predictor step satisfy for all 


Ui,j+j' = 0 (2.24) 

= 0 (2.25) 

= 0 (2.26) 

= 0 (2.27) 

5l+^'kk+r,, + h+r,j) = 0 (2.28) 


where K is defined in fl2.13l) . Then the finite volume scheme preserves the lake at rest 
and the jet in the rotating frame. 
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Proof. Since this argument is already standard for the lake at rest, we only sketch it 
briefly for the jet in the rotating frame. Let us study conservation of momentum in the 
^/-direction over cell Using fl2.24l) - fl2.26p . fl2.28l) . and the discrete product rule 

S{ab) = (5(a)p(6) -|- ii{a)5{b) (2.29) 

it is straightforward to show that the sum of the flux differences and the source term 
vanishes: 


- (/«.) 




Si’+’'Chuv) + Y‘^i- (Sp-’(hv- + W) + 9 ipHi + U)) . ( 2 , 30 ) 


Now = 0, and 




.y ^ , ^y y ^ J .^y 

= 0 . 


Dropping the corresponding terms in 02.301) . we obtain 

- (to)”- = Y sp-’s(h + b + U) 

i' 

= Y^‘' ■5“'V 

i' 

= 0 . 


This is the desired well-balanced property for the ^/-momentum. The x-momentum and 
the balance of mass can be treated analogously. ■ 


3 The well-balanced approximate evolution opera¬ 
tors 

The predictor step in the FVEG scheme will be based on exact and approximate integral 
representations of solutions to the linearized shallow water equations. We begin this 
section by formulating the exact integral representation. A few clarifying remarks should 
help the reader to understand the structure of this representation. Details of the derivation 
are given in Appendix We proceed to derive two approximate integral representations. 
For the hrst order scheme the approximate evolution operator for the piecewise 

constant data is used. For the second order method the continuous bilinear recovery 
R is applied. In the case of discontinuous solutions slopes in R are limited yielding a 
discontinuous piecewise bilinear recovery R, cf. Section 4 as well as [26]. The predicted 
solution at the quadrature nodes on the cell interfaces at the half timestep is obtained by 
a suitable combination of and 

E^nU" := + ES;?(1 - A4)ur 


( 3 . 1 ) 











Lukacova, Noelle, Kraft. 


Revised version, May 31, 2006 


where = l/4(f/j+ij+2f/jj+f/j_ij); an analogous notation is used for the y—direction. 
It has been shown in [2^ that the combination fl3.ip yields the best results with respect 
to accuracy as well as stability among other possible second order FVEG schemes. It 
is particularly important that the constant evolution term ~ corrects 

the conservativity of the bilinear recovery and hence of the intermediate solutions along 
cell-interfaces. If it is not used the scheme is second order formally, but unconditionally 
unstable, cf. the FVEG-B scheme |27j . 

Finally we show that approximate evolution operators lead to well-balanced two-step hnite 
volume schemes for the lake at rest and the jet in the rotational frame. 

3.1 Exact integral representation 

We believe that the most satisfying methods for evolutionary problems are based on 
the approximation of evolution operator or at least its dominant part. For the two- 
step FVEG method, we use two fundamental evolution operators. One of the steps is 
the classical hnite volume update for the cell averages and uses the integral form of the 
conservation law. Its well-balanced properties have been established in Theorem 12.II The 
other step, which precedes the hnite volume update, is needed to predict point values 
in order to evaluate huxes on cell interfaces. It is here that the classical bicharacteristic 
theory comes into play. It provides exact integral formulae for point values of solutions 
to multidimensional hyperbolic systems. 

Let P := {x,y,tn+i/ 2 ) be one of the quadrature points where the hnite volume huxes will 
be evaluated, and let w = {h, u, v) be a suitable local average of the solution around P. 
We will derive an exact integral representation of the solution of the linearized shallow 
water equation at P. Similarly as in fl 2 . 2 p . the linearized system in primitive variables 
reads 


Wt +Ai{w)w^ +A2{w)wy = t{w), (3.2) 

where the Jacobian matrices Ai and A 2 are dehned in fl2.3p . 

The homogeneous part of 03.21) yields a hyperbolic system. Fix a direction angle 6 with 
corresponding unit normal vector (cos 6, sin 6). The matrix pencil A = A{w) = cos 6Ai + 
sin 6 'A 2 , has three eigenvalues 


Ai = cos 0 M-|-sin 0 h — c, 
A 2 = cos^M-|-sin6*h, 

A 3 = cos 9 u + sin 9 V + c, 


and a full set of right eigenvectors 


-1 

ri = I gcos9/c 
g sin 6 */c 


r2 = 



rs = 


gcos9/c 
gsin9/c 


where c = 1 /^ denotes the wave celerity. The eigenvalues Ai^s correspond to fast waves, 
the so-called inertia-gravity waves, whereas slow modes are related to A 2 . Analogously 
to the gas dynamics the Fronde number Fr = \u\/c plays an important role in the 
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classification of shallow flows. The shallow flow is called supercritical, critical or subcritical 
for Fr > 1, Fr = 1, and Fr < 1, respectively. 

Applying the theory of bicharacteristics to the linearized system fl3.2p yields an exact 
integral representation of the solution. Since the computations are closely related to [27] , 
we summarize the key steps only briefly and refer to Appendix A for further details. 

• Fix a point P = {x,y,tn + t), t = 

For each spatial direction (cos( 6 *), sin( 6 *)), 6* G [0,27r) apply the corresponding one¬ 
dimensional characteristic decomposition to two-dimensional system fl3.2p . 

• Integrate the resulting equations along each bicharacteristic curve from time to 
time tn + T. 

• Integrate the resulting equations over all direction angles 9. This gives a represen¬ 
tation formula for the solution at the point P. 


P = {x,y,tn + t) 



Figure 1: Bicharacterestics cone. 


The EG integral representation derived in Appendix A then reads, cf. flA.40p - flA.42p 



(3.3) 


2 vr g 


/ 

Jtr 




I 



v{Q) cos 9 — u{Q) sin 9 j d6*dt. 
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1 g 


u{P) = 7 :u{Qo) + — / — (Q) cos6* + M (Q) cos^^ + n (Q) sin6*cos6*d6* 


27r 


hxiQo) + bx{Qo )) dt 


rtn~\~P 


J j cos^ 6* + sin6*cos6*) d6*dt 

1 1 

/ 


(3.4) 


r-27r 


tn tn + T — t Jq 

tn+T f Ptn+T [>2 tt 


u{Q) cos 29 + v{Q) sin 29 ) d6*dt 


/ 


^ i, X 

Z-n V Zji •J U 


v{ff) cos^ 9 — u{Q) sin 9 cos 9 ) d6*df. 


= o^(Qo) + 


1 9 


27r 


— ^h (Q) sin 9 + u (Q) sin 9 cos 9 + v (Q) sin"^ 9 d^ 


hyiQo) + by{Qo )) dX 


rtn+r 


— ^ J (^xiQ) S'ln 9 cos 9 + by {Q) siir^ 9 ] d 9 di 

ptn-\-T ^ 


(3.5) 


-27r 


27r ./f„ + T — t Jq 


f 


'tn 

rtn+T 


u{Qo)dt + 


271 


u{Q) sin 29 — v{Q) cos 29j d9dt 

tn+T 1‘2-K 

v{Q) sin 9 cos 9 — u{Q) sin^ 9 ) d6*dX 


Evolution takes place along the bicharacteristic cone, see Fig. 1, where P = {x,y,tn + t) 
is the peak of the bicharacteristic cone, Qq = {x — ur, y — vr, tn) denotes the center of the 
sonic circle at time tn, Qo = {x — u(tn + t — i),y — v(tn + r — i),i), Q = {x — u(tn + r — 
i) + c{tn + T — i) cos 9,y — v{tn + T — i) P c(tn + T — t) siu 9, i) stays for arbitrary point on 


the mantle and Q = Q{t) 

tn. • 


t — t-n 


denotes a point at the perimeter of the sonic circle at time 


At hrst view, it may seem difficult to interpret the terms in fl3.3p - fl3.5p . However, even if 
one is not familiar with bicharacteristic theory, there is a simple explanation of all terms, 
which we sketch in the following paragraph. 


з. 2 A simple interpretation of the EG integral representation 

Our interpretation of the EG integral representation is based on a comparison with the 
approximate representation of the solution which results from a Taylor expansion along 
the streamlines. 

The streamlines of the linearized shallow water equations are given by {{x(t),y(t),t)\x = 

и, y = h}. As in Section 2.2, let tp ■= tpt + itPx + v'Py be the material derivative of a 
function (p : Z) —)■ M. 








Well-balanced EG schemes. 


Revised version, May 31, 2006 


11 


Then the linearized shallow water equations (3.1) reduce to 

h = —h{ux + Vy) (3.6) 

ii = -Kx (3.7) 

V = —Ly. (3.8) 

This implies that 

h h{xix hy) hi^Kxx T ^yy) (^*^) 

U = -kx = gh{Ux + Vy)x - g{uhx + Vhy)x + fly (3.10) 

V = -Ly = gh{Ux + Vy)y - g{ubx + vby)y - fKx, (3.11) 


so the exact solution derived in fl3.3p - fl3.5l) can be approximated (up to 0{At^)) by 
fhi\ /h\ /h{Ux + Vy)\ f ^ h{Kxx +Lyy) \ 

I 1 • I ^ 1 ^1 j T 2 I gb-(^Ux T ^y^x gi^^bx T vby^x T f^y )• 

\^1 / /Qo \ /Qo \gh{Ux Vy)y — g(ubx vby)y — f Kx /q^ 

(3.12) 

Now we will indicate how to compare the RHS of the integral representation fl3.3l) - fl3.5p 
with that of (13.121) . We will show that they agree up to terms of 0{At‘^). For this we 
need the following Taylor expansions on the sonic circle and the bicharacteristic cone. For 
simplicity we introduce the notation 


a := cos^, /3 := sin 6*. 

Lemma 3.1. (Taylor expansions on the sonic circle and the bicharacteristic cone) Let 
(f : G™'+^(Z);M) and cpo := (p{Qo), (fo ■= g^iQo)- Then 


- 27 r 


k+l<m 1 r2n 


— / a^lS^ipde = — 

^ k,l=0 


kill 


(cr)'^+' + 0(Ar+^) (3.13) 


rtn+T r'iw 


2'kt 


a 


k +l<m 1 r2tr ol+n J a ^ r^n+r 

'I3’'^mdi= ^ J£Jo “ ^ ‘*''1 

k,l=0 


{c{tn + T-t)Y^^ d’).d[(podt 


kill 


T 


+ c>(Ar+^). 


Moreover, = Q if either k or I are odd integers, and 


1 

p27r 

1 


aMO = — 


JO 

27 r 

1 

p2-K 

1 


a^dO = — 


JO 

27 r 


r-27r 


kde = -, 
2 


- 27 r 


kde = -, 


c2tt 


2 ir L 8 


(3.14) 

(3.15) 

(3.16) 

(3.17) 


Proof, can be obtained by a direct evaluation. 


Using this Taylor expansion together with appropriate smoothness assumptions it is an 
elementary exercise to prove that 


j = |^tij(P) + 0(At='). 


(3.18) 
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3.3 Approximate evolution operators 

In the previous section, we could give the integral representation fl3.3p - fl3.5p a straight¬ 
forward interpretation by deriving it from a Taylor expansion along the streamlines. 

However, it would be far too expensive to evaluate fl3.3p - fl3.5p at each quadrature 
node of the two-step hnite volume scheme. In the present section we derive the crucial 
approximations of fl3.3p - fl3.5p . leading to an efficient and accurate algorithm. This 
approximation comes in two steps. First we derive in Lemma 3.2 a suitable approximation, 
which contains all terms necessary for the balance between the pressure terms and the 
sources, i.e. = 0 = Ly. This approximation is still continuous. Afterwards, we apply 
a special numerical quadrature to approximate the mantle integrals (i.e. time dependent 
integrals), in order to obtain approximate evolution operators which are explicit in time. 

For the present paper, we are interested in second order schemes. It is therefore sufficient 
that the predictor steps is hrst order accurate, i.e. accurate up to terms of order C>(Af^). 
In order to obtain a fully explicit hrst order approximation of {h,u,v){P), we would like 
to convert the mantle integrals on the LHS of fl3.3l) - fl3.5p into integrals over the sonic circle 
Sq. One possibility would be, analogously to |23], to use the simple rectangle rule 

-— / / (p{i,6)d6di= — / (p{tn,0)d6 + 0{At). (3.19) 

Jtn Jo Jo 


In the second step we can further eliminate derivatives over the sonic circle by means of 
the per-partes formula, cf. Lemma 2.1 in |23] 


CT 1 

Y Y 


-27r 


{(fix + Jjy)d9 


1 

— / {aip + l3ij)de + 0{At^). 
27r Jq 


(3.20) 


However, it has been shown in |23lE7ll2H] that the application of classical quadrature 
rules, such as the rectangle rule in fl3.19p . are not well suited for approximation of dis¬ 
continuous waves, which may propagate along the mantle of the bicharacteristic cone. It 
resulted in a reduced stability range of the FVEG. In particular, if the mantle integrals 
are approximated by the rectangle rule the CFL stability number was 0.63 and 0.56 for 
the hrst and second order FVEG scheme, respectively; it is the so-called FVEG3 scheme, 
cf. [21] . In the recent paper [27] new quadrature rules have been proposed for the mantle 
integrals. For example, if / = f{x) is a piecewise constant function, then it was shown in 
[27] that 


-27r 


27r 


f{Q )cos 9d9 + 




r2n 


- 27 r 


27r 


tn + T -t 


f{Q) COS 9di = 


27r 


f{Q) sgacos9d9, 

(3.21) 

an analogous relation holds for f{Q) sin 6*. Similarly, the quadrature rules for bilinear data 
have been derived, cf. Lemma A.l in the Appendix of 


As a result new approximate evolution operators evaluate exactly each planar wave prop¬ 
agating either in x— or y— directions and increase the stability range of the FVEG scheme 
substantially yielding the GFL number close to 1. 


Before we apply the quadrature rules proposed in [27], we approximate and simplify the 
exact integral equations fl3.3p - fl3.5p . This is done in Lemma 13^ Our strategy is to drop 
as many of the second order terms as possible, but to keep all those terms which enter 
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the balance of convective fluxes and source terms, i.e. = 0 = Ly^ and are therefore 
needed for well-balancing. The remaining terms will be reformulated or approximated up 
to the order 0{At‘^) in such a way, that the above mentioned quadrature rules from [27] 
can be applied. Thus, we keep the balance conditions for source terms and at the same 
time we will be able to approximate all resulting mantle integrals in a stable way. 

Lemma 3.2. The following operator is a first order approximation of the exact integral 
eguations 


- 27 r 


h{P) = -b{P) + ^l {h{Q) + h{Q))--{u{Q)cose + v{Q)sine)d0 

Jo 9 


T.l 


tn+T 


r-27r 


u{Q) cos 9 -|- v{Q) sin 9 ) dBdt 


tn + T -t Jo g 

tn+r /•2 tt 

ubx{Q) + vby{Q) ] dBdt + O , 


(3.22) 


r-27r 


U{P) = 


2ti 


— -K [Q) cos 9 + u{Q) cos^ 6* -I- n (Q) sin 6^ cos 9 d9 






r‘27r 


tn + r -t Jo c 


K{Q) cos 9 dBdt 


1 




1 


r-27r 


tn+T -t 


u{Q) COS 26* -I- v{Q) sin 26* j dBdt + O (At^) , 


(3.23) 


1 1 

V (P) = — / —L (Q) sin6*-|-M (Q) sin6*cos6*-|-n (Q) sin^ 6* d6* 

27r ./n c 


'0 


+ 0 "^ (Qo) “ 


I'tn+T 


o 2 -k 


2 

-/ 

2ir A 


27r 


tn + T -t Jo C 


L{Q) sin 6* dBdt 


(3.24) 


tn+T 


r2TT 


tn + T -t 


u{Q) sin 26* — v{Q) cos 26* j d9dt + O (At^) . 


The proof of this lemma is postponed to Appendix B. Now, in order to obtain time explicit 
approximate evolution operators we approximate time integrals in fl3.22p - fl3.24p . 

The only integral which is not of the form fl3.2ip is the last term in fl3.22p . Here we apply 
the rectangle rule fl3.19p and get 

j (ubxiQ) + vbyiQ)) d9di = — j (ubxiQ) + vbyiQ)) d9 + O {At^). 



Moreover, for the special case when the bottom topography slopes bx, by are approximated 
by a piecewise constant functions, which is the case of our bilinear recovery, for example, 
we can evaluate this term exactly. Note, that u = const., v = const., and b = b{x,y) does 
not change in time. Thus, we have 


1 



(^ubx{Q) + vby{Q)^ d9dt 


T 



{ubx{Q) + vby{Q)) d6*. 


All other integrals are of the form fl3.2ip . so we can apply numerical quadratures from 
f27\ . They will be used separately for constant and bilinear approximations. Thus, using 
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fl3.2ip we get, analogously to [27], the approximate evolution operator using the 

piecewise constant approximate functions 


r-27r r 


h(P) = 


u(P) 

v{P) 




{h (Q) + b{Q)) “ “ (Q) sgn(cos 6) + v (Q) sgn(sin 6)) 


de 


r-27r 


+ 7^- / {ub^{Q) + vby{Q)) d6 

27r Jq 

1 


r-27r 


27r 


27r 


r-27r 


— -K (Q) sgn(cos 9) + u (Q) ( cos"^ ^ d- -^] + v {Q) sin 6 cos 6 


— -L (Q) sgn(sin 9) + u (Q) sin 9 cos 9 + v (Q) ( sin 9 + - 

c \ ^ 


d9 


d9. 


(3.25) 


For the piecewise bilinear ansatz functions the justihcation of the mantle integrals ap¬ 
proximation is more involved. The reader is referred to [27] for more details. Applying 
the approximations from [27] the approximate evolution operator reads 


»2tt 


h{P) = -b{P) + h{Qo) + b{Qo) + -J^ {h{Q)-h{Qo)) + {b{Q)-b{Qo))d9 


2tt 


1 

TT 


0 


—u{Q) cos 9 -|- -v(Q) sin 9 
.9 9 


1 I 


r2TT 


d9 + — {ub^{Q) + vby{Q)) d9 

27r Jq 


u{P) = u{Qq) -/ -iF(Q) cos^d^ 

Jo c 

r2w 

2 


'0 


3u{Q) cos"^ 9 -|- 3v{Q) sin 9 cos 9 — u{Q) — -u{Qq) 


d9 


v{P) = v{Qo)- 


1 1 


TT 


L{Q) sin 9d9 


r-27r 


3u{Q) sin 9 cos 9 -|- 3v{Q) sin^ 9 — v{Q) - v{Qo) 


d9. (3.26) 


The approximate evolution operators fl3.25p . fl3.26p together with the finite volume update 
02.141) define the FVEG schemes. We will summarize the algorithm in the Section 4. 

Let us pause for a moment and discuss the possible advantages of using the approximate 
evolution operators instead of the Taylor expansion 03.121) in the predictor step entering 
the finite volume update 02.14p . 


• The above approximate evolution operators does not rely upon any derivative of 
the unknowns {h,u,v), while the Taylor expansion uses hx,hy,Ux + Vy. Therefore 
03.251) . 03.261) are less dependent upon the reconstruction of the piecewise constant 
solution than 03.12p . 

• The integrations along the cone in 03.3p - 03.5l) take into account the whole domain 
of dependence. This may result in a more robust algorithm, in particular for dis¬ 
continuous solutions. 
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3.4 Well-balanced property of the approximate evolution oper¬ 
ators 


The aim of this subsection is to verify that the approximate evolution operators fl3.25p . 
fl3.26p are well-balanced for the lake at rest as well as for the jet in the rotational frame. 
The proof of Theorem 13.11 below is based on the sufficient conditions for well-balancing 
formulated in Theorem 12.11 

Theorem 3.1. Suppose that the reconstructions at time satisfy for all {x,y) 


u^{x,y) = 0 (3.27) 

dyv^{x,y) = 0 (3.28) 

p.^v'^{x,y)dyh'"{x,y) = 0 (3.29) 

d^K'^{x,y) = 0 (3.30) 

dyL^{x,y) = 0. (3.31) 


Then the approximate EG predictor steps defined by fl3.25p . fl3.26p satisfy the conditions 
for well-balancing of Theorem 12.11 

Therefore, the FVEG schemes based on the above approximate evolution operators are 
well-balanced for the lake at rest and the jet in the rotational frame. 


Proof. We prove here that the approximate evolution operator for piecewise constant 
data see fl3.25p . satishes conditions fl2.24p - fl2.28p of Theorem 12.11 The proof for the 

approximate evolution operator for piecewise bilinear data 77“*"', see fl3.26p . is analogous. 


First we use conditions fl3.27p - fl3.3ip to simplify the approxmate evolution operator 
flT^ . Due to and flT^ 


r-27r 


r-27r 


— / M (Q) sgn(cos6*)d6* = — / n (Q) sgn(sin6*)d6'= 0 


27r 


27r 


Similarly 


r-27r 


r-27r 


r-27r 


— / 77 (Q) sgn(cos6*)d6* = — / L (Q) sgn(sin 6*)d6* = — / n (Q) sin^ cos6*d6* = 0, 


27r 
while 


0 


27r 


27r 


T Th (Q) (sin^» + 1) d» = + v% 


Due to ([S3ID, (1221 and ([221 


r-27r 


r-27r 


— / vby{Q)d9 = -— vhy{Q)d9 = 0. 


271 


277 


(3.32) 


Using the above identities in fl3.25p gives the simplihed approximate evolution operator, 
valid for the jet in the rotational frame 


ft (F) = -b(P) + ^ / (ft (Q) + i>(0))d« 

u (P) = 0 

viP) = !(«'■ +«*). 


(3.33) 
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From this, fl2.24p . fl2.25p and fl2.28p follow immediately. To verify fl2.26p . we set Pq := 
P — (0, 0, r) (the projection of P onto the plane t = and compute 

V iP) d,h (P) = f. (P) dy |-6(P) + ^ £’(h W) + 6(0))d«} 

= -i, (P) a,HPy) 

= -V (P) dy mp„) + b(p„)) - ;i(p„)) 

= f. (P) dyhiPy) 

= 0 . 


It remains to prove that fl2.27p holds. Since K = g{h + h — V), 

K{P) = g(h{P) + b{P)-V^+^/^{P)^ 

= ^ W) + KQ) - r”(Q)) + ^ r v’\Q)de - g V'”+‘''"(P), 

= T f K’'{Q)d0 + ^ I r”(Q)de - g V'“+'''2(P). 

Jo Jo 

Differentiating this equation with respect to x and applying fl2.12p implies 


- 27 r 


dggk{P) = ^1 d,K^{Q)d9 + ^ I d,V^{Q)d9-gdg,V^+^/\P) 


r-27r 


27r 

/ 


r-27r 


27r Jo 

g,R 

= f 


yR + yL 


2tt _ 
v^{Q)d9 - fv{P) 

R + v^_ 

J o 


which is the well-balanced condition (I2.27p and concludes the proof. 


4 Summary of the FVEG algorithm 

In this section we summarize the main steps of the FVEG method by presenting the al¬ 
gorithm for the first and second order scheme including the effects of bottom topography 
as well as the Coriolis forces. 


Algorithm 

1 Given are piecewise constant approximations at time h"-, n”-, i, j G Z, the 

bottom topography y), mesh and time steps h. At and constants g, /; compute 






n 

i'j 




9 


l'=lO 


rrn ^ ^ Kp-l + K' 

9 kri 


2 
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2 recovery step: 

If the scheme is second order, do the recovery step. For smooth parts of solntion 
apply the continuous bilinear recovery, cf. [26]. Possible overshoots on discontinuities 
are limited, e.g. by the minmod limiter; cf. [26]. This yields the piecewise bilinear 
approximations Rnu^, Rw", RnR^-, RnU'^-, RnV^- 

3 predictor step / approximate evolution: 

Compute the intermediate solutions at time level tn+ 1/2 on the cell interfaces by 
the approximate evolution operators. For the hrst order scheme use the approxi¬ 
mate evolution operator (I3.25p : the second order scheme is computed using 

both approximate evolution operators fl3.25l) as well as fl3.26p . cf. fl3.ip . 

Integration along the cell interfaces is realized numerically by the Simpson rule. 

4 corrector step / FV-update: 

Compute the Coriolis forces and the bottom topography at the intermediate time 
level tn+ 1/2 and at each integration points on cell interfaces, i.e. at vertices and 
midpoints: 


jn+ 1/2 


T 7-71+1/2 

'^7 + 1 / 2 ,^ 


r 7-71+1/2 
^A3j + l/2 


b{xk,y£), k = i,i±l/2,i 


f 




71 + 1/2 , 71 + 1/2 

C-1/2,^ + ^i+l/2,£ 


9 


71 + 1/2 , 71 + 1/2 

'^fcj-1/2 + '^kj+1/2 




= J,J±l/2; 


k = i,i ±1/2. 


Do the FV-update fl2.14p using the well-balanced approximation of the source terms 

dMU). 


5 Numerical experiments 

One interesting steady state, which should be correctly resolved by a well-balanced scheme, 
is the stationary steady state, i.e. h + b = const, and u = 0 = n. In this section we demon¬ 
strate well-balanced behaviour of the proposed FVEG schemes through several benchmark 
problems for stationary and quasi-stationary states, i.e. h + b ^ const, and u ^ 0 ~ n; 
see [HI m\ for related results in literature. Further, we present results for steady jets 
including effects of the Coriolis forces and show that the FVEG scheme is well-balanced 
also for this nontrivial steady state. At the end of this section we compare accuracy and 
computational time of the well-balanced FVEG method and the well-balanced second and 
fourth order FVM of Audusse et ah [2] and Noelle et ah [31]. 


5.1 One-dimensional stationary and quasi-stationary states 

In this experiment we have tested the preservation of a stationary steady state as well as 
the approximation of small perturbations of this steady state. The bottom topography 
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consists of one hump 


b{x) 

and the initial data are u{x, 0) 

h{x, 0) 


0.25(cos(107r(a; - 0.5)) + 1) if |a: - 0.5| < 0.1 
0 otherwise 


1 — b{x) + e if 0.1 < a; < 0.2 
1 — b{x) otherwise. 


The parameter e is chosen to be 0, 0.2 or 0.01. The computational domain is the in¬ 
terval [0,1] and absorbing boundary conditions have been implemented by extrapolating 
all variables. The gravitational constant g was set to 1 analogously as in [niEO]. It 
should be pointed out that the one-dimensional problems are actually computed by a 
two-dimensional code by imposing zero tangential velocity n = 0. 


Firstly we test the ability of the FVEG scheme to preserve the stationary steady state, 
i.e. the lake at rest case, by taking e = 0. In Table 1 the L^-errors for different times 
computed with the hrst order FVEG method, cf. fl3.25p . and with the second order FVEG 
method, cf. fl3.26p . are presented. Although we have used a rather coarse mesh consisting 
of 20 X 20 mesh cells, it can be seen clearly that the FVEG scheme balances up to the 
machine accuracy also for long time computations. 


Table 1: The L^-error of the well-balance FVEG scheme using 20 x 20 mesh cells. 


Method 

t = 0.2 

t = 1 

t = 10 

first order FVEG 
second order FVEG 

0 

1.67 X 10-1^ 

0 

1.11 X 10-1^ 

2.22 X 10-^“ 

4.27 X 10-1® 


Top surface at time t=0.0 Top surface at time t=0.7 




Figure 2: Propagation of small perturbations, e = 0.2. 

In Figure 2 the typical propagation of small height perturbations is shown at time t = 0.7. 
The solution is computed on a mesh with 100 x 5 cells and the height of the initial pertur¬ 
bation was e = 0.2. The initial disturbance generates two waves, the left-going wave runs 
out of the computational domain, and the right-going wave passes the bottom elevation 
obstacle. It is known that if the perturbations are relatively large in comparison to the 
discretization error a “naive” approximation of the source term, i.e. not well-balanced 
scheme, e.g. the fractional step method, can still yield reasonable approximations. How¬ 
ever, for small perturbations, i.e. e of order of the discretization errors, such a scheme 
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Fignre 3: Propagation of small perturbations, magnified view; e = 0.2 (left) and £ = 0.01 
(right). 

would yield strong oscillations over the bottom hump and the wave of interest will be lost 
in the noise, see [20] . 

In Figure 3 we compare results for water depth h at time t = 0.7 obtained by the hrst and 
second order FVEG methods using the minmod limiter and the monotonized centered 
limiter (denoted as MNG), respectively. In the left picture e = 0.2, the right picture 
shows results for e: = 0.01. The reference solutions was obtained by the second order 
FVEG method with the minmod limiter on a mesh with 10000 cells. For the hrst order 
scheme and the second order scheme with minmod limiter we can notice correct resolution 
of small perturbations of the stationary steady state even if the perturbation is of the 
order of the truncation error. The MNG limiter resolves the peak much more sharply, but 
overcompress the left-going wave. This is a well-known feature of compressive limiters, 
see e.g. the discussion in [331 1^ flS]. 

5.2 Two-dimensional quasi-stationary problem 

The second example is a two-dimensional analogue of the previous one. The bottom 
topography is given by the function 

0.8 exp (-5 {x - 0.9)^ - 50 (?/ - 0.5)^) (5.34) 


_ f l — h{x^y)+e if 0.05 < X < 0.15 
( 1 — b{x, y) otherwise 

= 'i;(x,|/,0) = 0. (5.35) 

The parameter £ is set to 0 and 0.01. The computational domain is [0, 2] x [0,1] and the 
absorbing extrapolation boundary condition are used. 

First, we take e = 0 and test the preservation of a two-dimensional lake at rest on a mesh 
with 20 X 20 mesh cells, see Table 2. Analogously to the one-dimensional case this steady 
state is preserved up to the machine accuracy. In the Figure 4 we present two solutions 
of a perturbed problem, e = 0.01, which are computed on a 200 x 100 grid (left) and on a 


h{x,y) 

and the initial data are 

h(x,|/,0) 

m(x,|/,0) 
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Table 2: The L^-error of the well-balance FVEG scheme using 20 x 20 mesh cells. 


Method 

t = 0.2 

t = 1 

t = 10 

first order FVEG 
second order FVEG 

2.35 X 10-1^ 

4.97 X 10-1^ 

5.09 X 10-1^ 

6.74 X 10-1^ 

5.02 X 10-1^ 

1.53 X 10-1® 


600 X 300 grid (right) by the second order FVEG scheme with the minmod limiter. Notice 
that the FVEG method correctly approximates small perturbed waves, the perturbation 
propagates over the bottom hump without any oscillations. Note that the wave speed is 
slower over the hump, which leads to a distortion of the initially planar perturbation. The 
perturbed wave runs out of the computational domain and the flat surface is obtained at 
the end. Our results are in a good agreement with other results presented in literature, 
cf., e.g., 

5.3 Steady jet in the rotational frame 

This is a classical Rossby adjustment of an unbalanced jet in an open domain, see e.g. [5]. 
The initial data are a rest state superimposed by a one-dimensioal jet, 

h{x, y, 0) = 1.0, u{x, y, 0) = 0, v{x, y, 0) = 2Nl{x), 

where the shape of the velocity v is given by a smooth profile 

, . {1 + tanh{4x/L + 2)) {1 — tanh{4x/L — 2)) 

^ {l + tanh{2)y 

with L = 2. We have used flat bottom topography b{x) = 0, the parameter of the Goriolis 
forces / and the gravitational acceleration g are set to 1. The nondimensional parameter 
representing the effects of Goriolis forces, the Rossby number Ro = = 1 and 

the Burgers number reflecting the nonlinear effects is Bu = = 0.25. The initial 

jet adjusts a momentum unbalance, which emits the waves, the so-called gravity waves, 
propagating out from the jet. The formation of shocks can be noticed within the jet core 
approximately at vr//, which is a half of a natural time scale Tf = 2n/f, see Figure 5, 
6 . As time is evolved the solution tends to the equilibrium state fv = gh^, which is 
a geostrophic balance as demonstrated in Figure 7. We can notice that even for long 
time simulations there are still small oscillations around the geostrophic equilibrium. As 
pointed out by Bouchut et ah [5] some wave modes with the frequencies close to / remain 
for a longer time in the core of the jet. Their analysis for a linearized situation shows that 
they correspond to the gravity wave modes having almost zero group velocity, and thus 
are almost not propagating. For another extensive study of the stability of jets, which 
gives interesting eigenfunctions similar to those in Figure 5 we refer to m- 

5.4 Accuracy and performance 

In this experiment we compare accuracy and computational time of the well-balanced 
FVEG, the second order well-balanced FV method of Audusse et ah [2] and its fourth 
order extension due to Noelle et al. [31]. We choose a fully two-dimensional experiment 
analogous to that of Xing and Shu [M] , but include moreover the effect of Goriolis forces by 
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Top surface at time t=0.6 


Top surface at time t=0.6 




Top surface at time t=0.9 


Top surface at time t=0.9 



Top surface at time t=1.2 


Top surface at time t=1.2 




Top surface at time t=1.5 


Top surface at time t=1.5 




Top surface at time t=1.8 



Top surface at time t=1.8 



Figure 4: Two dimensional quasi-st at ionary problem fl5.34p . fl5.35p . 


setting / = 10. The gravitational constant was set to g = 9.812. The bottom topography 
and the initial data are given as follows 

b{x,y) = sin( 27 ra;)-f cos( 27 r|/), 
h{x,y,0) = 10-|- exp(sin( 27 ra;)) cos( 27 rj/), 
hu{x,y,0) = sin(cos( 27 ra;)) sin( 27 r 2 /), 
hv{x,y,0) = cos( 27 ra;) cos(sin( 27 r?/)). 

The computational domain [0,1] x [0,1] was consecutively divided into 25, 50,..., 800 
mesh cells in each direction. We have compared solutions obtained by the second order 
FVEG scheme as well as by the second order and fourth order well-balanced FVM at time 
T = 0.05. For the second order well-balanced FVM of Audusse et ah the second order 
Runge-Kutta method was used for time integration, the third order Gaussian quadrature 
was used for cell-interface integrals of fluxes and the second order WENO recovery was 
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applied. The reference solution was obtained by the fourth order well-balanced FV method 
of Noelle et al. [31] . 

Tables 3 and 4 contain the errors and experimental order of convergence (EOC) for the 
FVEG for both CFL numbers 0.8 as well as 0.5, respectively. The well-balanced higher 
order directional splitting FVM is in general stable up to CFL=0.5. The errors for 
its second order version are presented in Table 5. We can indeed see that both methods 
are second order accurate in all components. Moreover, the second order FVEG scheme 
is almost 10 times more accurate than the second order FVM, see Table 5 as well as the 
left picture of the Figure 8. 

Figure 8 illustrates the GPU/accuracy behaviour graphically. We use the logarithmic 
scale on x—,y— axis. On the y— axis the errors in Erst component h is depicted. 
Errors in other components yield analogous results. On the left of Figure 8 a comparison 
between second order FVEG and FV methods are presented, whereas on the right we 
show the comparison between the fourth order well-balanced FVM of Noelle [21] and 
the second order FVEG scheme. The FVEG schemes yields on coarse meshes still more 
accurate solutions. In fact, for meshes up to approximately 100 x 100 cells, which are 
actually often used for practical computations, it is more efficient to use the second order 
FVEG scheme than the fourth order FVM. The superiority of the fourth order scheme is 
showing up on hue grids, see the right graph of Figure 8. 

We should point out that no attempt has been made in order to optimize the codes with 
respect to their GPU performance. Our extensive numerical treatment indicates that 
both well-balanced second order methods, the FVEG as well as the FVM are actually 
comparable with respect to their computational time. 


Table 3: FVEG scheme: Gonvergence in the L} norm, GFL=0.8 


N 

error in h 

EGG 

error in hu 

EGG 

error in hv 

EGG 

25 

1.04e-02 


3.56e-02 


8.52e-02 


50 

2.42e-03 

2.10 

8.71e-03 

2.03 

2.15e-02 

1.99 

100 

6.01e-04 

2.01 

2.23e-03 

1.96 

5.50e-03 

1.96 

200 

1.54e-04 

1.96 

5.76e-04 

1.95 

1.44e-03 

1.93 

400 

3.97e-05 

1.96 

1.47e-04 

1.97 

3.69e-04 

1.96 

800 

1.02e-05 

1.97 

3.71e-05 

1.98 

9.40e-05 

1.97 


Table 4: FVEG scheme: Gonvergence in the norm, GFL=0.5 


N 

error in h 

EGG 

error in hu 

EGG 

error in hv 

EGG 

25 

1.37e-02 


6.19e-02 


1.18e-01 


50 

2.80e-03 

2.29 

1.05e-02 

2.56 

2.33e-02 

2.34 

100 

5.23e-04 

2.42 

1.80e-03 

2.54 

4.25e-03 

2.45 

200 

1.04e-04 

2.33 

3.63e-04 

2.31 

8.12e-04 

2.39 

400 

2.45e-05 

2.09 

8.79e-05 

2.05 

1.80e-04 

2.17 

800 

6.14e-06 

1.99 

2.20e-05 

2.00 

4.36e-05 

2.04 
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Table 5: FV scheme: Convergence in the norm, CFL=0.5 


N 


FOG 

error in hu 

EGG 


FOG 

25 

4.53e-02 


2.13e-01 


3.40e-01 


50 

1.32e-02 

1.77 

5.57e-02 

1.94 

9.51e-02 

1.84 

100 

3.50e-03 

1.92 

1.42e-02 

1.97 

2.52e-02 

1.92 

200 

8.95e-04 

1.97 

3.58e-03 

1.99 

6.46e-03 

1.96 

400 

2.26e-04 

1.99 

8.96e-04 

2.00 

1.63e-03 

1.99 

800 

5.67e-05 

1.99 

2.24e-04 

2.00 

4.10e-04 

1.99 


6 Conclusions 


In the present paper we have developed a new well-balanced scheme within the frame¬ 
work of the hnite volume evolution Galerkin (FVEG) scheme. The scheme is applied 
for the shallow water equations with source terms modelling the bottom topography and 
the Coriolis forces. The key ingredient of this FVEG scheme is a new well-balanced ap¬ 
proximate representation of the solution, cf. Lemma 3.2, which together with a recent 
quadrature rule from 1271 leads to the multidimensional approximate evolution operators 
fl3.25p . fl3.26p . These approximate evolution operators are used in a predictor step. In fact 
we are predicting the solution at cell interfaces and do not need to use the hydrostatic 
reconstruction as it is done by Audusse at el. [2]. 

In the following correction step, which is the hnite volume update step, the source term is 
approximated in the interface-based way. We have proved that the lake at rest, the steady 
jet in the rotational frame as well as their combinations are preserved, cf. Theorems 2.1 
and 3.1. Numerical experiments in one and two space dimensions demonstrate correct 
resolution of these equilibrium states and of their small perturbations. For smooth solu¬ 
tions the accuracy of the well-balanced FVEG scheme is superior to that of a recent FV 
scheme while the GPU time is comparable. 

In future work we want to extend our well-balanced schemes to shallow water equations 
with nonlinear friction, which appears in oceanographic as well as river how modelling. 
Another challenge is presented by multi-layere shallow water models, which are important 
in oceanology, meteorology and climatology. 
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A Derivation of the exact integral equations 


Applying theory of bicharacteristics to fl3.2p we can derive exact integral eqnations in 
an analogous way as in |27]. In order to keep the presentation self-contained we briefly 
rewrite main steps of the derivation. 

Let R be the matrix of right eigenvectors corresponding to direction n := (cos( 6 '), sin( 6 *)). 
Its inverse reads 

/ —1 - cos 9 - sin 6 * \ 

1 / ^ ^ \ 

R = - 0 2 sin 6 * —2 cos 6 

^ \ 1 - cos 9 - sin 9 I 

Let us define the vector of characteristic variables v by 

V := R~^w. 


Multiplying system (13.2p by ^ from the left yields the following characteristic system 

dv ^ ~ 

- + = r, 

where 


Bi = 


Bo = 


u — c cos 9 

— 1 h sin 6* 

0 \ 

—g sin 9 

u 

gsin9 

0 

1 h sin 6^ 

u + c cos 9 / 

V — c sin 9 

1 h cos 9 

0 \ 

g cos 9 

V 

1 

o 

o 

0 

— 1 hcos9 

V + c sin 9 / 

ri \ 

1 

\-g{{-9hx + fv) cos9 - {gby + fu) sin^) 

To = ^’ 

\n)t = 

{—gbx + fv) sin 9 + {gby -|- fu) cos 9 

rs / 

\ 

l^{{-gbx + fv) cos9 - {gby + fu) sin^) 


r(n) = 


and the characteristic variables v are 


(A.36) 


v{n) = 



|(—h -I- ^ucos9 + |nsin 6 () 
= R^^(n)u= I M sin 0 — n cos 6 * 

\{h+ % cos 9 + sin 9) 


The quasi-diagonalised system of the linearized shallow water equations has the following 
form 


dv 

'm 


( u — c cos 9 

0 

0 ^ 

1 dx 1 

( V — c sin 9 

0 

0 \ 

0 

u 

0 

0 

V 

0 

V 0 

0 

■u - 1 - c cos 9 ) 

^ 0 

0 

h - 1 - c sin 6 * / 


dv 

— = s + r 
dy 

(A.37) 


s = 



/ 


\h{su).9^+ cos9^) \ 

g sin 9(^ -^)-g cos 9(^ - fa) 

^ k ax / J \ an oy ' 


dx dx ' ^ ^ dy 

fh{cos9^ - sm9^) 


with 
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Let us denote by Xi the Gth bicharacteristic corresponding to the Gth equation of system 
(IA.37p . The bicharacteristic xi is dehned in the following way 


dxe{s) f b\f \ 
ds 


where bj^, 6 ^^ are the diagonal entries of the matrices Bi, B 2 , respectively. The bicharac¬ 
teristics Xi create the surface of the so-called bicharacteristic cone, see Fig. 1, with the 
apex P = {x, y, tn + r) and the footpoints 


Qi{0) = (x — {u — ccos 9)T,y — {v — csin 9)T,tn), 

Q 2 = Qo = {x - UT,y -VT,tn), 

Qsid) = {x — {it + c cos 9)T,y — {v + c sin 9)T,tn)- 

Remember that r = At/2 in our case. Integrating each equation of flA.37H along the 
corresponding bicharacteristic from the apex P down to the footpoints we get 

ptn+T 

ve{P) = ve{Qe)+ se{Qe{i)) + ri{Qi{t))dt, £=1,2,3. (A.38) 

J tn 

Now multiplying flA.38p with R from the left and averaging over all directions we go back 
to the original variables w 


w{P) = 

1 72^ 
2vr Jo 


-vi{Qi{9), 9) V3{Q3{9), 9) 

lcos9vi{Qi{9),9) + sin9v2{Q2{9),9) + f cos 0 ^ 3 (( 53 ( 6 '), 6 ')) 
|sin 6 'ni(Qi( 6 '), 6 ') - cos9v2{Q2{9),9) + ^sin9v3{Q3{9),9)) 


d9 


1 r 2 n f -s'i{9)-r[{9) + s^{9)+r'^{9) \ 

+ W / ^ cos 9{s[{9)+r[{9)) + sin 9 {s 2{9)+r'2{9)) + ^ cos 9{s'^{9)+r'^{9)) d0, 

^ y ^ sin 9 {s[{9) + r[{9)) — cos 9 {s' 2{9) + r2{9) + ^ sin 9{s^{9) + r'^{9)) ) 

(A.39) 

where s^( 6 *) = si{xi{i,9),i,9)di is an integral along the £-th bicharacteric and the 

analogous notation holds for source terms r^. It should be noted that the source term 
in flA.3811 arrives from to the multidimensionality of the system, whereas the source term 
re is a physical source term. 

Now, we have Ai = —A 3 , Qi{9 -|- vr) = Q3{9) and the characteristic variables ve are 27r- 
periodic. Applying the Gauss integration, cf. (I3.20p . in order to avoid the derivatives of 
dependent variables appearing in s we can, after analogous computations as in [2311^ . 
reformulate the exact integral equations flA.39D in the following way 


h{P) = 


17^’" c c 

— / h{Q) - u{Q) cos 9 -n(Q)sin6'd6' 

27r Jo 9 9 

— — / -= / - (m(Q) cos6'-|-n(Q) sin^ ) d6*df 

27r tn P T — t Jo 9 

1 rtn+T |•2^T 

+—C I J \bx{Q) cos 9 + by{Q) sin 9) d9dt 

rtn+T /•27r 

\'^(Q) COS 6 * — m(Q) sin 6*) d^df, 

27r 9 Jtn Jo 


(A.40) 


27r 
1 cf 
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r-27r 


9 


(Q) cos 6 + u (Q) cos^ 6 + v (Q) sin 6 cos 9 dO 


- I 

2 / 

’27r 

1 




hx{Qo) + bx{Qo)) dt 

ctn+T I'2'K 


(A.41) 


t-n 


^tn-\-T 


bx {Q) cos^ 9 + by (Q) sin 9 cos 9 ) d6*dt 

1 




tn + r-t Jq 

v{Qo)di + 7 ^ 

Z7T 


u{Q) cos 29 + v{Q) sin 29 ) d6*dt 


rtn+T /•2w 


v{Q) cos^ 9 — u{Q) sin 9 cos 9 ) d6*df. 


v{P) = ^^(<50) + ^ 


l^2n 


— ^h {Q) sin 9 + u {Q) sin 9 cos 9 + v {Q) sin^ 9 d9 


- I 

2 / 

’27r 

1 

■il 


rtn+r 


hyiQo) + by{Qo) ) dt 

I'tn+T I'2'K 


(A.42) 


bx (Q) sin 9 cos 9 + by (Q) sin"^ 9 ) d6*dt 


r'^n~\~T 


r‘27r 


tn + T — t 




u{Qo)dt + 


271 


u{Q) sin 29 — v{Q) cos 29j d9dt 

tn+T 1‘2-K 

v{Q) sin 9 cos 9 — u{Q) sin^ 9 ) d6*dt 


Recall that Qq = {x — uit^ + r — t),y — v{tn + T — Q = {x — u(tn + r — t) + c(tn + t — 
i) cos 9,y — v{tn + T — i) + c{tn + T — i) sin 9, i) stays for an arbitrary point on the mantle 

and Q = Q{t) _ denotes a point at the perimeter of the sonic circle at time 

t — tn 


B 


Proof of Lemma 




Proof. We show here that the approximate integral equations fl3.22l) - fl3.24|) are consistent 
with the exact integral equations fl3.3l) - fl3.5p . i.e. flA.40p - flA.42p . In fl3.3p the integral with 
bottom topography terms can be rewritten using the polar-type transformation along the 
mantle of the bicharacteristic cone, i.e. xq = a:-|-r(cos6* — |), y^ = y + r{sm9 — |), where 
r = c{tn + r — i) is the circle radius at the time level i G [t„, + r]. Thus, we have 

^(r,0) = bx{xQ,yQ)cos9 + hy{xQ,yQ)sm9 - ^ 2/q) + |/q)) . 

Therefore, 


c 


c 


p 2 , 7 r 

In Jo 

0 /•27r 


bx{Q) cos 9 + by{Q) sin^ j d^df 


(B.43) 



CT Jo 




ftn+T /‘2 it 


ubx{Q) + vby{Q)d9di 


'0 


"" d / 1 
dr V 27r 

2w 


c2tt 


_(_/ bd9]dr + —J^ 


"tn+T 1‘2'K 


I'tn+T 


= 77- / b{Q)d9 - b{P) + — 


271 


>0 


271 


ubx{Q) + vby{Q)d9dt 


ubxiQ) + vby{Q)d9dt, 


2it 



















Well-balanced EG schemes. 


Revised version, May 31, 2006 


27 


which yields the corresponding terms in fl3.22p . 

Further, we show that the integrals in fl3.3p containing the Coriolis forces are of order 
0{At‘^)] note that r = At/2. Applying the rectangle rule in time and the Taylor expansion 
from Lemma 13.11 in the center of the sonic circle Qq yields 


rtn+T /'2 tt 


r2n 


v{Q) cos9d9dt = T / v{Q)cos9d9 


r-27r 


= r / {v{Qo) cos9 + CTVx{Qo) cos^ 9 + CTVy{Qo) cos9sm9 + 0{At))d9 

Jo 

= 0{At^) (B.44) 

with an analogous approximation for the Coriolis forces in ?/-direction. Together with 
(IB.43P and flB.44p this yields the approximate integral equation fl3.22p . 

In the equation fl3.4p for velocity u we apply for the mantle integrals containing the bottom 
elevation terms the rectangle rule in time and the Taylor expansion over the center Qq of 
the sonic circle S'o at time which lead to 

tn-\-T ‘2'K 

9^ r 


1 


g 



bx{Q) cos9 -|- by{Q) sin^ ) cos9 d6'dt = + 0{At‘^). (B.45) 


To complete we eliminate the derivative by replacing the term bx{Qo) by its average over 
the sonic circle Sq and applying the Gauss theorem 


c2tt 


bxi.Qo') 


bx{Q) dxdy 0{At^) = 


'•S'o 


TlC^At^ 

which after substitution into flB.45p yields 


TTcr 


b{Q)cos9d9 + 0{At^), (B.46) 


1 


tn+T 2ir 


9 



bx{Q) cos 9 -|- by{Q) sin6') cos6' d^dt = 


C 271 


r-27r 


b{Q) COS 9 d9 + 0{At^ 


tn 0 


(B.47) 

Rewriting the Coriolis forces terms using their primitives we obtain analogously to flB.45p 
and flB.46p 

rtn+T /•2 tt 

v{Q) cos^ 9 — u{Q) sin 9 cos 9 ) d^dt (B.48) 


271 


tn 

271 
9 1 


rtn+T p2Tr 


r-27r 


Vx{Q) cos 9 — Uy{Q) sin6') cos 9 d6* dt 


= -— / V{Q)cos9d9 + 0{Af). 

c27i Jo 

This balances together with flB.47p the analogous term with h{Q) cos 9 in fl3.4p . The 
integral along the middle bicharacteristic 


1 

2 


rtn+T 


hxiQo) + bx(Qo) — —v{Qo)\ dt — ^ 
9 J ^ 


rtn+T 


hxiQo) + bxiQo) — VxiQo)) dt 


can be approximated in a similar way as flB.46p applying the Gauss theorem at each inter¬ 
mediate circular section at t along the mantle of the bicharacteristic cone. Substituting 
into fl3.4p gives fl3.23p . Approximation fl3.24p for the velocity v is obtained in an analogous 
way as fl3.23p . ■ 
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height; t=0 


height: t=1.257 




Figure 5: One-dimensional Rossby adjustment problem, time evolution of water height. 
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height; t=0 


height; t=1.257 



height; t=2.5152 


height: t=3.7692 



Figure 6: Rossby adjustment problem, time evolution of water height, two-dimensional 
graphs. 
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1=136.9017 1=138.2325 




1=140.0023 



Figure 7: One-dimensional Rossby adjustment problem at different times, geostrophic 
balance. 


FVEG (2nd order, CFL=0.5), FVM (2nd order, CFL=0.5) 


FVEG (2nd order, CFL=0.8), FVM (4th order, CFL=0.5) 




log(CPU Time) 


log(CPU Time) 


Figure 8: Efficiency test: error over the CPU time; the second order FVEG and second 

order FV schemes (left) as well as the fourth order FVM and the second order FVEG 
scheme (right). 










































